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Abstract We discuss capability of Smooth Particle Hydrodynamics to represent ade- 
quately the dynamics of self-gravitating systems, in particular for what regards the quality 
of approximation of force fields in the motion equations. When cubic spline kernels are 
used, we find that a good estimate of the pressure field cannot be obtained in non uniform 
situations using the commonly adopted scheme of adapting the kernel sizes to include a 
fixed number of neighbours. We find that a fixed number of neighbours gives the best 
approximation of just the intensity of the force field, while the determination of the direc- 
tion of the force requires a number of neighbours which strongly depends on the particle 
position. A good balance between quality of the approximation and computer time con- 
sumption is obtained assuming the latter about 2 times larger than the first one. We 
suggest a procedure which is massively parallelizable. With a suitable choice of the kernel 
sizes the amount of computations required is less than twice the one required by the com- 
mon SPH scheme at the same level of approximation, so our recipe (when parallelized) is 
convenient. 

1. Introduction 

The subject of this paper is a numerical technique which is largely used in 
fluid-dynamical simulations: the so called Smooth Particle Hydrodynamics 
(SPH). 

This particle method was firstly introduced by Lucy [1] and Gingold and 
Monaghan [2] to simulate nonaxisymmetric fluids in Astrophysics, and it 
became popular due to its robustness and easy implementation. 

Basically, the idea of the method is to consider the fluid as an ensem- 
ble of (smooth) particles, anyone being representative of a piece of fluid. In 
this scheme, each particle is mathematically characterized by a "kernel" (a 
symmetric, regular non-negative function centered on the particle position) 
which carries information on the average values of dynamical and thermody- 
namical quantities, as well as on their gradients. The size of the kernel may 
depend on the local density of the particles. 

Each particle moves in the force field generated by the whole particle 
system, while the associated quantities evolve under their suitably regularized 
laws (a more detailed description of this method can be found in [3]). 

It is well known that necessary requirements for a good numerical method 
are compatibility with "analytical" equations, convergence to the right solu- 
tions and stability. At present, these properties have not yet been completely 



2 



proved for SPH. Few results are available; for instance Oelschlager [4] and Di 
Lisio [5] have obtained convergence results for a particular equation of state 
(EOS), P oc p 2 , in the free and self-gravitating cases, respectively. 

A more general convergence result, for a freely moving polytropic com- 
pressible fluid, can be found in Di Lisio, Grenier, Pulvirenti [6], [7]. 

These results seem to show that the SPH method has the desired char- 
acteristics of convergence, at least for polytropic EOS. 

Of course, one of the most important requirement for the method to work 
well is providing a good evaluation of the force field. 

Usually, in the simulations, the regularization kernels have compact sup- 
port (to reduce the computational cost) and their widths are such that every 
particle interacts with an almost constant and a priori fixed number of neigh- 
bours. It is found that this evaluation strongly depends on the size of the 
kernels. 

In this paper we attempt to give a better way to establish, through a 
quantitative check of the error, a correct size of the kernels for a good evalu- 
ation of the force field, in the class of the cubic spline kernels. Actually, the 
spline kernels introduced and discussed by Monaghan and Lattanzio [8], are 
the most used in numerical simulations. 

In what follows we shall denote with V or V r the usual simbolic vector 
of components J^, J^, while the symbol * represent the usual convolution 
operator in H 3 The points in the space will be denoted by r and r = |r|. 

2. SPH approximation of the fluid equations 

2. 1 The free case 

A particle method to simulate the dynamical and thermodynamical evo- 
lution of a fluid corresponds to a direct discretization of the Lagrangian set 
of fluid-dynamical equations. In the case, for instance, of a barotropic com- 
pressible fluid without body forces this set writes as 

f f + P V " v = ° 

< f + ¥ = ° (!) 
P = P(P)- 

In the SPH scheme the fluid, whose total mass is M, is represented by an 
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ensemble of N particles with masses mi (i = 1, . . . , N) such that J2i m i — 
M, centered at points ri; {ri} is a set of space points sampling the volume 
occupied by the fluid. To each particle is associated a density function called 
"kernel" . The most commonly used kernel is the third order spline function 

(see [8]) 

1 ( 1 - (3v 2 /2) + (3t» 3 /4) if < v < 1 
0(r,/i) = — (2-vf/A if l<v<2 (2) 

[ elsewhere 

where v = r/h. The width of the kernel, 2h, is the radius of the spherical 
"fluid particle" ; this radius may change particle by particle. 

In SPH the pressure field per unit mass in the Euler equation is usually 
approximated by 

(f) (r ( ) = E(^ + ^W ( - rj ,fc,, (3) 

\ P J SPH j=l \ Pi Pj / 

where 

N 

Pi = p(ri) = m j(p( r i - r i, hi) (4) 

3=1 

is the SPH approximation to the fluid mass density. 

The density (|j) evolves in time because the particles move under the 
Newton law 

* = -(—) (rt) i = l,---,N (5) 

V P ) SPH 

which is the SPH version of the momentum equation in the system ([[]). 



2.2 The self-gravitating case 



SPH methods developed mainly in Astrophysics, due to its ability to 
handle fully 3D situations. Most of astrophysics simulations should neces- 
sarily include self-gravity. The system of equations underlying the motion of 
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barotropic self-gravitating fluids is a little more complicated than (1): 

I + PV-V = 
— 4- — - VU = 

dt p (6) 
V • VU = -4nGp { ' 

P = P{p) 

where G is the gravity constant. Due to the spherical simmetry of the kernels 
and Newton theorems, the SPH approximation of the gravitational force field 
is simply: 

(V£/) 5 PH ( r *) = " E ^g \ Yi ~ V3 \ /' ri ~ rjl <f>{x, K)x 2 dx (7) 
~{ \Ti-rj\ 6 Jo 

This SPH approximation is quite accurate, and numerical simulations show 
that the gravitational field is always fitted better than the pressure one (be- 
cause of that the sum in (^) is usually substituted by a further approxima- 
tion as in the P3MSPH code suggested by Monaghan and Lattanzio [8] or in 
TREESPH code proposed by Hernquist and Katz [9]). It is easy to convince 
oneself about that when considering the different nature of the pressure and 
gravitational fields. The first depends on the local distribution of the matter, 
and so it is strongly subject to the local fluctuations of the particle distribu- 
tion, while the second rather depends on the global matter density. This is 
the reason why, in a particle scheme, it is easier to fit the gravitational force 
field than the pressure one. 

So, in this paper we will focalize our attention on the SPH approximation 
of the pressure field, only. 



2.3 Evaluating the error 



The h value determines the number of particles which give an effective 
contribution in the sums (|3]) and (^), that is the number of neighbouring 
particles that effectively interact (touch) with the given one. This means 
that the effective pressure force strongly depends on h: a small h means a 
steep (f) profile in ([|) and a small number of particle contributing to the sums, 
and viceversa for large h. It is easy to verify the two limits: 

(r) = lim (— - J (r) = 0. (8) 

SPH h -*°°\ P J SPH 



lim f — — ] 

h-oV p ) 
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In the common case of a polytropic EOS (P oc p 7 ) the second limit holds 
just for 7 > 2/3. It is so quite clear how the kernel size h is crucial to 
give a good SPH approximation of the pressure force field so that the par- 
ticles, fluid-representatives, move on the correct trajectories (which are the 
characteristics of the fluid). 

A quite natural way to determine a "best" local value of h is through the 
minimization of the absolute value of the relative error 



E r {v- h) 







W 




(¥) 


M 





(9) 



At a given particle configuration, for fixed r, E r depends on h only. 

If we set the unitary vectors n e (r) = (— )/| — | (the exact direction) and 
(^j-)sph/\^j-\sph (the approximated SPH direction) the relative 



error 



can be rewritten as 



E r {Y- h) 



n e - n a + 



1 — 1 


— 1 — 1 


1 p 1 


1 p 1 



SPH 



| VP I 
I p I 



-n. 



(10) 



This obvious identity suggests us to split the relative error into a direction 
error E dir and a modulus error E mod . Respectively 



E dir (r; h) = |n e - n | , -E mo d(r; h) 



1 — 1 — 


1 — 1 


1 p 1 


1 p 1 



SPH 



VP I 



11^ 



3. A model for comparing exact and SPH pressure fields 



To test the sensitivity of the SPH approximation of the force field on the 
kernel width, we will consider particle distributions sampled from the family 
of Plummer's density laws [10]: 



p(r) 



1 + (r/r c y 



5/2 • 



(12) 



The law (|12D is often used in theoretical Astrophysics because it is a good 2- 
parameter (po and r c ) fit to many self-gravitating objects (globular clusters, 
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roundish elliptical galaxies, etc.). Here it is used just as a reference density 
law apt to simulate a range of different density and pressure gradients. To do 
this one can change the central density po an d the lenght scale r c ; increasing 
r c means dilute the distribution (smoothing the gradients). 

We fix the total mass of the sphere, M so that the central density p is 
linked to the characteristic length scale r c which is our second free parameter. 
As we said, in the SPH scheme the density is represented by an ensemble of 
N particles with masses m 8 and positions r$ (i = 1, . . . , N). In this scheme 
the pressure field in the position can be approximated by (f|). The particle 
distribution sampling the regular density ( |i2| ) is obtained choosing first a set 
of random particle positions according to a uniform probability law in a 
sphere and then replacing every Fj by rj such that 

J%(s)s 2 ds = T f. (13) 

In our simulations we shall use 10, 000 particles of individual mass = 
M/10,000. 

An interesting quantity associated to h is the number of particles which 
give an effective contribution in the sums (|3|) and (^), that is the number of 
particles that effectively interact with the given one. In practical simulations, 
this number is usually set to a constant value, the same for every particle. 

Due to the spherical simmetry of our model, the directional part of the 
relative error Edir is a decreasing function of h, and is obviously minimized 
taking h as large as possible (which is computationally expensive), while the 
error in the absolute value E mod is not a monotonic function of h (see eq.(|)). 
So, it is quite natural to study the modulus part of the error separately. 



3.1 The results 

We have studied the minimization of the relative error E ri as a function 
of h (assumed to be the same for each particle) for three different values of r c 
in fljjp, namely 0.25, 0.5, 1. For 7 the values 7/5, 5/3 and 2, of clear physical 
interest, are considered. 

For each case the minimization of E r as function of h has been performed 
on a set of 100 different particle positions. To do that we evaluate the 
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pressure field on each particle using the same kernel for all the particles, 
then we modify the size of this kernel to obtain the best result. 

Fig.s 1 and 2 show (for some r c and 7) the minimum relative error E r as 
function of the distance from the centre and the number of particles (neigh- 
bours) within the kernel of any of the 100 particles, respectively. 

These figures seem to show that the best results (low value for E r coupled 
with a small number of neighbours) have been obtained in the middle radius 
zone. The minimum error E r is, in the average, a slightly decreasing function 
of r. 

To study better the behaviour of the E r function, we have done a sepa- 
rate statistics of the modulus (E mo d) and directional (E^ir) errors. We have 
considered 30 particles in three different regions, around r/r c = 0, 1, 2, as- 
suming, in the state law, 7 = 5/3. To give a visual idea of the quality of 
the SPH approximation, in Fig.s 3, 4, 5 the SPH pressure field is compared 
with the exact one. For graphical convenience, for each particle we have 
plotted the exact vector and a vector whose modulus is given by the best 
SPH approximation (that given by E r minimization) at an angle, clockwise 
oriented, given by the scalar product between the exact, n e , and SPH, n a , 
unitary vectors in space (cp. Sect. 2.3). The exact vector points to the 
centre; in the figures the centre is assumed to be at an infinite distance. It is 
quite evident the better field direction approximation for larger r c and r/r c 
(i.e. less steep gradients). 

To give a more precise measure of the errors, in Fig.s 6, 7, 8 we have 
plotted the relative error E r and the modulus error E mod . The number of 
neighbours and, to give a measure of the direction error E^ir, the cosine of 
the angle G between the SPH and exact vector fields (cos = n e • n a ) have 
also been plotted. Fig.s 6 — 8 show that the field is well approximated in 
the middle zone while a very good approximation of just the direction of 
the pressure field is obtained in the outer one. The best balancing of the 
two errors have been obtained in the middle zone. In the inner zone, where 
the pressure gradient should almost vanishing, the approximation is critically 
dependent on sampling problems (see discussion in Sect. 4.). 

The different radial behaviours of the errors E r , E mo d and E^ir, suggest 
us to split the evaluation of SPH approximation in two parts: 
first we fix h to evaluate at best the modulus of the pressure field, then we 
enlarge it to obtain an acceptable direction of the field with a computationally 
reasonable number of neighbours. 
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To this scope, for several particles positioned at different radii, we have 
plotted the errors E r , E mo d and Edir as functions of h. These plots, Fig. 9, 
show that the best error in the pressure field absolute value (E mod ) is reached 
with a smaller h value than when the global error (E r ) was minimized. 

Taking into account these results, we suggest to estimate separately the 
modulus and the direction of the pressure field using two different kernel 
widths, as aforementioned. The first width, say h mo d, is obtained through 
the minimization of the modulus error E mo d, the second, h d i r , is taken pro- 
portional to the first in a way to keep the error in direction Edir below an 
acceptable limit. 

To set such proportionality constant, we evaluate, for any particle in the 
three regions considered before, the quantity Edir as function of the variable 
h/hopt, where h opt is the width obtained minimizing the modulus error E mod . 

The plots in Fig. 10 have been obtained averaging the E dir (h/h opt ) func- 
tions found for each particle. The lenght of the vertical straight lines is 
proportional to the dispersion around the mean value. 

This figure shows that it is sufficient to double the h opt value to get a good 
(within 10%) approximation of the direction of the field. So a reasonable 
choice of the proportionality constant is 2. 

In Fig.s 11, 12, 13 comparisons of the results obtained following this 
recipe compared with those obtained by the simple minimization of E r are 
shown. In these figures the number of neighbours corresponding to h opt 
and the percentage global error E r are plotted for all particles in the three 
zones before considered for r c = 0.25, 0.5, 1, respectively, and 7 = 5/3. 
The figures show that h opt is associated to an almost constant number of 
neighbour particles, which corresponds to about 10% of the total number. 
This behaviour occurs in all the cases here considered. 

Note that usual SPH implementations evaluate the pressure field using 
kernels with a size such that the number of neighbouring particles is set to a 
given constant. Thus, we tested the computational cost of the implementa- 
tion here proposed by comparing the CPU times (on a DEC Alpha 200 4/233) 
required to evaluate the pressure field by the usual and our SPH method. 

We have fixed the number of neighbouring particles to compute the ab- 
solute value of the pressure gradient to 10% of the total, as suggested by our 
E mo d minimization. This fraction corresponds to a certain kernel size, and to 
evaluate the direction of the field we double this size. In the usual method 
just one kernel size is used (here we choose the one giving 10% neighbours). 
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As axpected, the computational cost of a double-kernel evaluation of the 
pressure field is ~ 1 + 2 3 times the single-kernel one, which is the price of the 
good improvement of the error. 

For the sake of a good balance between computational weight and pre- 
cision we tried hdi T = \^h mo d, i.e. the size needed to almost double the 
number of neighbours. This case requires a computational cost which is of 
the order of 1 + 2 times the single-kernel method with h = h mo d, with a lower 
error level; to have a similar approximation we verified that h mo d is allowed 
to be < h thus to give a ratio of CPU times less than 2. Being this part of 
SPH evaluation wholly parallelizable explains why in this way we can obtain 
a better approximation waiting the same solar time. 

4. Comparing SPH and Monte-Carlo methods 

The basic idea of the SPH method is to "reconstruct" a function /, 
whose values are known in some "particle" points, by interpolation (in fluid- 
dynamical applications the particles are distributed to fit the fluid mass den- 
sity). The way SPH does that is through a Monte-Carlo approximation of the 
convolution integral f *<p, where is a regular kernel function. The function 
/ * is close to / if is close to the Dirac's delta. In the same way it is 
possible to reconstruct V/. Indeed, by noting that (V/) * <fi = f * (V0), we 
can evaluate the latter integral, instead of the first, through the knowledge of 
/ only. A similar procedure may be used for any derivative of / with respect 
to space variables. 

Monte-Carlo methods are the most studied and the most advanced parti- 
cle methods in numerical analysis, so it is quite natural to compare the SPH 
approximation method with direct Monte-Carlo ones. 

While we refer the reader to a forthcoming paper [11] for a more detailed 
and accurate analysis of the link between SPH and Monte-Carlo approxima- 
tion methods, for the purposes of this paper we limit here ourselves to show 
the difference between the SPH and the Monte-Carlo interpolations for the 
pressure force field only. 

As we said in previous Sections, the pressure force field VP/p is the 
most delicate term to handle in fluid-dynamical simulations. In SPH, VP/p 
is usually approximated by eq.(|). The Monte- Carlo rule, choosing p as 
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distribution function of the particles, gives 



\ P ) MC N i P ^) 

where we assume hi = h for any i. The average value, over the particle 
configuration space, of the sum in the right hand size of eq. ([14]) gives the exact 
pressure field (VP/p)(r). On the contrary, averaging the SPH approximation 
sum ([|) we find 



■ , , (r) = — (r) + / l^pr - fg±) Vp(r')0(r - r', h)dh', 

\ P J SPH \ P J J \P W P( r )> 



/VP\ \ /x /VP\ , s rfP(r) P(r') 

-J Vp(_r' )(p[r — r', h) 

(15) 

the value of the integral being zero if P oc p 2 , otherwise it goes to zero with 
h whenever the functions VP/ p 2 and Vp are regular enough (for example if 
they are differentiable) . This means that the SPH scheme (without artifi- 
cial viscosity) may not be apt to represent extreme situations, as when, for 
example, shock fronts occur (i.e. when a discontinuity on the density occurs). 

Now, it is easy to show that the difference A between the Monte-Carlo 
and SPH sums writes as 



1 E 



+ 



| E Vp(r,)0(r - r 4> h) - i £ ^V r 0(r - r„ h) (16) 

and its average value is not zero, as we said. 

This term is too difficult to handle in general situations. But we can 
study it in some simple cases. For example, for a fluid which state law is 
P oc p 2 , then eq. ([l6|) reduces to 

A = |E^ryV,(p(r^(r-r 4 )). (17) 

In a situation not far from spatial uniformity (| Vp|/p << 1), eq.(|TT|) reduces 
to: 

A ^EV ri 0(r-r 4 ), (18) 
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whose averaged value is the vector — Vp. On the other hand, if we suppose 
|Vp|/p» 1 

A 4?lr*-"> < 19 > 

holds, and the averaged value of the right hand term is the vector Vp. 

The first case should be the better one, anyway when the density profile 
is quasi-constant, the fluctuations due to the particle fit amplify the relative 
error. 

Let us consider now the state law P oc p 2+a where a > 0. When 
ln 2 (p)max(l, p~ a ) << 1 we have P ~ p 2 (l + aln(p)), and (|T^|) writes as 

A -|E^ryV, (p(r l )0(r-r i )) + 
^ £ V, pn(p(r 4 ))0(r - *)] + ^E^" *)V r , ln 2 (p(r 4 )). (20) 

i i 

Let us study again the two extreme cases |Vp|/p << 1 and |Vp|/p >> 1. It is 
easy to convince ourselves that the previous considerations still hold, but only 
if | lnp| is not too great. This means, for example, that SPH approximation 
is not expected to be accurate near the boundaries of the fluid. 

In conclusion, the best SPH fit (i.e. comparable with Monte-Carlo one) 
for the pressure field can be obtained when both the following conditions 
hold: 

i) the density is not too large or too small and ii) it has a profile along a 
well defined direction with small gradients. 

Moreover, we notice that the previous analysis does not depend on the 
overall mass distribution but only on its local behaviour. Then the conditions 
i) and ii) are valid criteria in very general situations. 

Let us now apply previous considerations to the spherical model studied 
in this paper. The results of our simulations, shown in Sect. 3, confirm 
the expectations based on the previous analisys: the best results on the 
evaluation of the pressure field are obtained in those zones of the model 
where the conditions i) and ii) approximately hold, i.e. in the central ones 
as suggested by the density law (|P2"D and 



— *T^A < 21 > 

p (r 2 + r~) 
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Moreover, comparing models with different characteristic radii we found that 
the precision increases as r c grows, as suggested by 

max^ocl. (22) 
p r c 



5. Conclusions 

In this paper we have studied the error in the SPH evaluation of the pres- 
sure field for the set of polytropic compressible fluids in spherical simmetry. 
The spherical simmetry allows us to compare SPH evaluation with a simple 
model, saving the confidence that the main results are still valid in more 
general and realistic situations, as we have shown with analytical consider- 
ations in Sect. 4. First we determined the kernel sizes which minimize the 
global error E r (see eq.(|])); they depend on the particle positions becom- 
ing narrower where the density is lower, due to that in external regions the 
right (radial) direction of the field is easily obtained. This latter is linked to 
that a good approximation of the pressure field is harder to be obtained in 
very central and peripheral regions. Actually, in very central zones, where 
|VP|/p ~ and |Vp|/p « 1, to get an acceptable fit of the pressure field 
direction so large kernel sizes are nedeed to make the SPH pressure gradient 
underestimated. On the contrary, in the outer regions, the radial direction 
of the field is well fitted even with small kernel sizes (i.e. few neighbours): 
this implies a certain overestimate of the pressure gradient. 

All this has suggested us to study separately the behaviours of the mod- 
ulus and directional error, as defined by eq. (|lT|). 

Given that the larger the kernel the better the approximated direction, we 
minimized the error in modulus E moi i rather than the global one. 

The obtained kernel width slightly depends on the particle position thus 
an almost constant number of neighbours is needed to perform the best 
approximation. In our simulations this number is found to be 10% of the 
total independently of the particle positions and the characteristic radius of 
the model. 

Moreover, a good (giving an error less than 10%) direction of the pressure 
field is obtained doubling that kernel size. Unfortunately, this evaluation 
requires a CPU time almost 8 times the one required by the single-kernel 
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evaluation. We have obtained a good balance between the quality of approx- 
imation and the CPU weight fixing the 'direction' kernel to be \/2 times the 
'modulus' one. 

So, a good prescription to estimate the pressure field in SPH is to evaluate 
the absolute value of the field with a kernel size giving a constant number of 
neighbouring particles, while the direction is estimated using another kernel 
size as described before. 

The CPU time required by this procedure is less than two times that 
required by the usual (single-kernel) SPH evaluation of the pressure field. 
This recipe becomes convenient when a parallelized code is used. 
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Figure captions 



Fig. 1 

For the Plummer's sphere with r c = 0.25, and polytropic exponent 7 = 7/5 (solid 
line) and 7 = 2 (dashed line) panel (a) shows the minimum relative error E r (in percent) 
vs. the radial coordinate and panel (b) the corresponding relative percentage number of 
neighbouring particles N n (see text). 

Fig. 2 

As in Fig.l but r c = 1. 

Fig. 3 

For the polytropic exponent 7 = 5/3, the exact (vertical arrows) and approximate 
VP/p fields are shown in three different radial zones, for r c = 0.25 . 

Fig. 4 

As in Fig. 3 but for r c = 0.5 

Fig. 5 

As in Fig. 3 but for r c = 1. 

Fig. 6 

For the polytropic exponent 7 = 5/3 and r c = 0.25, we show, from bottom up: 

- the minimum relative error E r in percentage; 

- the corresponding modulus error E mo( i in percentage; 

- the cosine of the angle O between the unitary vectors n a and n e (see text); 

- the number of neighbours N n in percentage to the total; 

for all the 30 particles in each of the three chosen shells. Solid line refers to the inner 
(r/r c = 0) shell, dashed line to the middle (r/r c = 1) and dot-dashed to the outer 
(r/r c = 2). 

Fig. 7 

As in Fig. 6 but for r c — 0.5 

Fig. 8 

As in Fig. 6 but for r c = 1. 
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Fig. 9 

For the polytropic exponent 7 = 5/3 and characteristic radius r c = 0.5, the values of 
the errors E r (solid lines), E mod (dashed lines) and E dir (dot-dashed lines) are plotted 
as functions of h/r c for 10 particles positioned at different distances from the centre (as 
labelled). 

Fig. 10 

For the polytropic exponent 7 = 5/3 and r c = 0.25, 0.5, 1 the averaged (over all the 
30 particles in each shell) directional error Edir as function of h/h op t is plotted (see text). 
The lenghts of the vertical straight lines are proportional to the dispersions of the data 
around the mean values. 

Fig. 11 

For the polytropic exponent 7 = 5/3 and r c — 0.25 the percentage number of neigh- 
bours N n and the percentage total error E r are shown for all the particles in the shells 
when our double-kernel method is applied (solid lines) and when E r is minimized (dashed 
lines) for all the 30 particles in the shells (as lebelled in the right panels). In each case the 
horizontal straight line correspond to the average value over all the 30 particles. 

Fig. 12 

As in Fig. 11 but for r c = 0.5 

Fig. 13 

As Fig. 11 but for r c = 1. 



17 



Toward the centre 




A 



Toward the centre 



nf f f f f if tinff fill hiUhUii: 




r/r = 1 

/ c 



A I I 



A '< » 



| I I 
I 

I 



A A 



A A ♦ 
A 

A 



A /I 



ii A h ' l ^ 



r/r =0 

/ c 



Toward the centre 



!! I H H 1 1 U i U I it i 1 1 I H I II it I I 



4 i J I 4 111 I A 



r/r =2 

/ c 

IMIMUHUIUl 



r/r =1 

' c 



1 1 

4 t 



4 * 



It 



It 



I 4 



M M M 



III! 



i \ 



ii ii 



r/r =0 

/ c 




10 20 30 



particle's label 




10 20 30 



particle's label 




10 20 30 



particle's label 



100 
50 


100 
50 


100 
50 


100 
50 


100 
50 





i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

h/r h/r 



150 



•I— I 

T3 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 




1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



r/r c =l 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



r/r=2 




-i i i i i i i i i i i i i i i i i i i 



1234 1234 1234 



h/h 



opt 




1234 1234 1234 



h/h 



opt 



150 i_i i i i i i i i i i i i i i i i i i i 1 1 1 i i i i i i i i i i i i i i i i i i 1 1 1 i i i i i i i i i i i i i i i i i i 



100 
50 





: r =0.25 

c 

-I I I I I I I I I I I I I I I I I I I 



I I I I I I I I I I I I I I I I I I I 




-I I I I I I I I I I I I I I I I I I I 



1234 1234 1234 



h/h 



opt 




10 20 30 10 20 30 

particle's label particle's label 



